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SUPPRESSION  OF  RIVER  ICE 
BY  THERMAL  EFFLUENTS 

George  D.  Ashton 


PART  I.  UNSTEADY  SUPPRESSION  OF  RIVER 
ICE  BY  FULLY  MIXED  THERMAL  EFFLUENTS 

Introduction 

Imposition  of  a  dam  and  a  reservoir  on  a  river  gen¬ 
erally  results  in  the  release  of  water  during  winter 
periods  at  a  temperature  above  that  which  would  or¬ 
dinarily  be  experienced  without  the  reservoir.  As  a 
result,  the  ice  cover  is  suppressed  in  the  reach  down¬ 
stream  of  the  dam.  The  present  analysis  examines  this 
ice  suppression  in  an  unsteady  sense  by  including 
both  the  effects  of  unsteady  variation  of  the  meteorolog¬ 
ical  variables  and  the  storage  and  release  of  the  energy 
associated  with  the  melting  and  freezing  of  the  ice 
cover.  The  effluent  from  the  reservoir  is  assumed  to 
be  fullv  mixed  ovei  the  initial  cross  section.  In  Part 
II  of  this  report,  the  effects  of  a  side  discharge  of 
thermal  effluent  will  be  considered  so  as  to  include  the 
effects  of  lateral  mixing. 

Previous  analytical  work  on  the  problem  is  sparse. 
Dingman  et  al.  (1968)  analyzed  the  steady -state  case, 
that  is,  assumed  constant  meteorological  conditions 
and  chose  as  the  criterion  for  the  location  of  the  ice 
edge  downstream  the  point  where  the  water  temperature 
was  0°C.  They  well  recognized  the  limitations  of  their 
steady-state  model  (Weeks  and  Dingman  1972).  Paily 
et  al.  (1974)  considered  a  similar  problem  but  with  a 
step  increase  in  the  effluent  temperature  and  inclusion 
of  the  effect  of  longitudinal  dispersion.  Their  criterion 
for  ice  edge  location  was  also  the  location  of  the  O^C 
isotherm  for  open  surface  conditions.  Examination  of 
their  results  also  shows  that  inclusion  of  the  longitu¬ 
dinal  dispersion  term  has  small  effect  compared  with 
simple  nondispersive  routing. 

In  contrast  to  these  analytical  studies  Donchenko 
(1978)  pointed  out  that  the  ice  edge  responded 


strongly  to  changes  in  the  meteorological  variables  and 
varying  discharges,  in  one  case  fluctuating  over  a  range 
of  80  km.  Donchenko  also  pointed  out  the  necessity 
of  considering  the  hydrodynamic  stability  of  the  leading 
edge  in  locating  the  position  of  the  ice  front. 

Governing  equations 

For  a  flow  with  temperature  uniformly  mixed  over 
the  depth,  with  no  transverse  velocity  variations,  and 
steady  flow,  the  governing  partial  differential  equation 
is 


(D 

where  D  is  the  depth  of  flow,  7"w  is  the  water  tem¬ 
perature.  U  is  the  mean  velocity  in  the  x  (longitudinal) 
direction,  p  is  the  density,  Cp  is  the  specific  heat  ca¬ 
pacity,  z  is  the  transverse  coordinate,  and  E x  and  £, 
are  respectively  the  longitudinal  and  transverse  mixing 
coefficients.  <}>  is  the  heat  flux  at  the  top  surface  and 
there  is  assumed  to  be  no  flux  at  the  bottom.  If  we 
assume  D,  U,  p  and  Cp  are  constant,  then  the  equation 
may  be  written  in  the  form 
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We  now  neglect  longitudinal  mixing  and  obtain 

OLa.  W iL*  JL  le,  — _  (3) 
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It  we  further  assume  that  is  fulls  mixed  across  the 
width,  then 
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I  malls,  it  a  lagrangian  viewpoint  is  adopted  we  mas 


Outline  of  analysis 

A  definition  sketch  is  presented  in  I  igure  I. 

When  the  top  surface  of  the  flow  is  open  to  the 
atmosphere,  o  where  is  the  heat  flux  front 
the  watei  surface  to  the  atmosphere  above.  is  a 
complicated  function  which  has  a  number  ol  components 
depending  upon  ait  temperature,  wind  velocity,  humidity , 
cloud  cover,  time  of  day  and  other  factors.  I  cm  present 
purposes,  we  will  assume  the  wind  speed  l  „  and  ait 
temperature  l  t  are  dominating  and  express  bs  the 
simple  relationship 

(rw-r*> 

where  is  a  heat  transfer  coefficient  that  is  a  function 
of  the  factors  noted  above.  As  a  practical  matter  we 
will  use,  in  the  numerical  simulation  which  follows,  a 
relationship  of  the  form 


wheie  a  has  dimensions  of  W  m'*'  C  and  b  has  dimen¬ 
sions  of  W  m' 's'Y  1 .  Depending  upon  the  asailabilits 
ol  input  data,  one  mas  alwass  calculate  an  equivalent 
value  ol  a,  b  oi  /rw4  bs  enetgy  budget  analsses.  In 
general,  a,  b  and  are  functions  of  time  only ,  since 
the  lengths  osei  which  the  analysis  is  cairied  suit  are 
generally  much  shorter  than  corresponding  spatial 
variations  ol  the  energy  budget  sanables.  However,  out 
main  purpose  here  is  not  to  examine  letinements  in 
energy  budget  calculations  but  rather  to  examine  the 
implications  id  the  neglect  of  the  unsteadiness  of  the 
total  flux  0  and  the  enetgy  required  to  melt  (oi  fiee/e) 
an  ice  cover . 

When  the  tlow  is  ice-covered,  b  but  now  the 
heat  flux  from  the  water  to  the  ice  depends  upon  the 
flow  variables.  In  particular,  we  take  the  heal  flux 
from  the  water  to  the  ice.  yS^,  as 


wheie  /i„,  is  a  heat  tianstei  coefficient  and  /  is  the 
melting  point  (/,„  0"C).  Hosed  conduit  turbulent 

heat  tianstei  correlations  are  geneiallv  of  the  form 


n  TTi  U  ) 


where  K  is  the  hvdtaulic  radius  (m),  (  M  is  the  average 
flow  velocity  (ms'1),  pw  is  the  water  density  (kgm'M, 

H  is  the  dynamic  viscosity  (kg  m  s  '),  ('  is  the  specific 
heat  (/  kg  1  V' ),  and  is  the  thermal  conductivity 
of  the  water  (W  m  ’  ('  is  an  empirical  coefficient 

i«n  the  order  of  0.023  (see,  e.g„  Kohsenow  and  Choi 
I  dpt)  for  smooth  sui  faces.  When  the  watei  is  above 


J 


freezing  and  the  flow  is  turbulent,  relief  features  (ice 
ripples)  form  on  the  underside  of  the  ice  cover  which 
increase  the  value  of  C  by  up  to  50%  (Ashton  and 
Kennedy  1972).  Evaluating  the  properties  at  0°C  and 
taking  R  =  £>/ 2, 


where  Cwi  =  1622  W  s0-8  m'2'6  °C'  for  C  =  0.023  and 
Cwi  =  2433  W  s0-8  m'2-6  0c’  for  C  =  0.0345. 

Equation  5  may  now  be  integrated  to  yield,  in  co¬ 
ordinates  moving  with  the  flow, 
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[  pCpD  J 

_  KMo)] 
7*.  0-Tm  ~  [  pC,D  J 


open  surface  (11) 


ice  covered  ( 1 2) 


where  the  initial  condition  is  taken  as  Tw  =  Tw  0  at 
t  =  t0  and  x  =  xQ.  Relative  to  coordinates  fixed  in 
space  such  that  dx  =  Udt  the  corresponding  equations 
are 


open  surface  (13) 


rw-ra=exp[ 
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PCPUD  J 

Tw~Tm  = 

pwi  (^oH 

w,  0“'  m 
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ice  covered.  (14) 


Equation  14  is  shown  in  Figure  2  for  typical  parameter 
values. 

The  melting  and  thickening  of  the  ice  cover  is  governed 
by  the  energy  balance  at  the  water/ice  interface 

0i-0wi=.piX^  (15) 

dt 

where  0(  is  the  heat  flux  by  conduction  through  the  ice 
cover  (W  m'2),  pt  is  the  density  of  ice  (kg  m'3),  X  is  the 
heat  of  fusion  ( /  kg'1 ),  tj  is  the  ice  thickness  (m),  and 
dr\tdt  is  the  rate  of  thickening.  The  conductive  flux 
01  through  the  ice  cover  is  treated  in  a  quasi-steady  man¬ 
ner  by  assuming  a  linear  temperature  profile  through 
the  ice  thickness.  Thus, 

.  _*|(rm-7s)  nr;* 


where  kf  is  the  thermal  conductivity  of  the  ice  (W  m' 
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X-X„  (km) 

Figure  2.  Downstream  attenuation  of  water  tempera¬ 
ture  beneath  an  ice  cover. 


°C'1 ),  Ts  is  the  top  surface  temperature  of  the  ice  cover, 
and  we  require  that  Ts  <  0°C  because  of  the  state  re¬ 
lationship.  In  the  absence  of  evaporation  or  conden¬ 
sation  on  the  top  surface,  0j  =  0jj  where  is  the  heat 
flux  from  the  surface  to  the  atmosphere.  0U  may  be 
calculated  in  a  manner  similar  to  0WJ  through  intro¬ 
duction  of  a  heat  transfer  coefficient  h b  applied  to  the 
difference  Ts-Ta  in  the  form 

07> 

Using  0 j  =  0U  allows  us  to  eliminate  Tf  between  cq  16 
and  1 7  and  results  in 


0i  = 


7m -Tt 

JL  +  J. 


(18) 


We  could  also  in  a  similar  manner  add  the  effects  of  a 
snow  cover  by  introducing  an  additional  term  rj i/ki  in 
the  denominator  of  cq  1 8. 

Substitution  of  eq  18  in  cq  15  then  results  in 

d9) 

JU_L  dt 

*i  Ttjj 

which,  if  Ta  and  Tw  are  constant  in  time,  may  be  readily 
integrated.  They  seldom  are,  however,  so  in  the  numer¬ 
ical  analysis  eq  19  will  be  integrated  numerically  in  step¬ 
wise  fashion  in  the  form 
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Location  of  ice  edge 

Before  examining  the  effect  of  vary  mg  7"  ,  it  is  of 
interest  to  examine  the  criterion  for  the  location  of  the 
ice  edge  undei  conditions  of  constant  Tt.  Previous  in¬ 
vestigators  have  genet  alls  uses)  the  location  of  the  0"l 
isotherm  predicted  t>v  analysis  01  the  heat  loss  (torn  the 
open  surface  (essentially  the  analvsis  leading  to  eg  1  3) 
with  various  modifications  from  the  present  analvsis  to 
consider  different  formulations  for  or  to  include 
the  effects  vvl  longitudinal  dispersion,  f  ront  eg  1  3  h\ 
the  present  analvsis. 
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f  iqurr  J.  Comparison  of  distune -es  do*  mtneum  to  ciiqe 
of  ice  coier  tor  the  criteria  ot  equilibrium  heut  tluses 
jrtd  0  C isotherm  locution .  us  u  function  o!  the  esufer 
remperuture  uir  temperuture  rutio. 


Since  the  heat  balance  at  the  surface  determines 
whether  or  not  ice  will  form  at  that  surface,  it  is  a  bet¬ 
ter  indication  of  the  location  of  the  ice  edge  than  is 
the  location  of  the  O'C  isotherm.  I  he  condition  is  ob¬ 
tained  bs  assuming  incipient  ice  formation  corresponds  to 
dr?  Jt  0  when  t?  =  0.  f  rom  eg  1 11  this  is  at  a  tem¬ 
perature  7,  t  given  h\ 

(22) 

Substituting  into  eg  1 3  the  corresponding  location  is 
given  bv 


(>-.Vn)l  -oCAV 


ft  is  now  necessary  to  assume  that  h  which 

is  reasonable  at  the  ice  edge  since  the  top  surface  of  ihe 
ice  evfge  may  be  considered  wet.  The  ratio  of  the  dis¬ 
tances  bv  the  two  criteria  is  then  given  b\ 


7j,  v  0  C.  Ihe  ratio  is  presented  in  I  igurc  3  as  a  func¬ 
tion  of  7*  q/(-7j)  for  various  values  ot  7tu  h^t.  The 
ratio  is  small  corresponding  to  small  salues  ot  7*  0  -7. 
(either  large  initial  water  temperatures  ot  small  sub- 
frec/ingair  temps'tatures)  and  large  salues  of  hu 
(laige  wind  speeds  ot  slow  flow  velocities).  Negative 
values  ot  the  ratio  correspond,  ot  course,  to  cases  where 
v  'vhich  physically  mas  lx-  interpreted  as  those 
conditions  tor  which  ice  ssill  thicken  became  the  heat 
loss  to  the  atmosphere  is  greater  than  the  heal  delivered 
to  the  undersurface  bv  the  (low.  to  put  the  situation 
into  better  perspective,  I  igurc  4  presents  salues  ot  the 
ratio  for  particular  values  of  7W  O'^wa-^ivi  and  74. 

Which  of  the  two  criteria  for  the  ice  edge  location 
is  operative  depends  on  the  seguence  of  aii  tempeiatuies. 
Oearls ,  water  will  not  free/e  until  its  temperature  has 
decreased  to  O'C.  Thus,  unless  a  prior  ice  cover  is 
present,  the  most  upstream  location  of  the  ice  edge  is 
given  bv  the  O'C  isotherm,  corresponding  to  the  criterion 
represented  bv  eg  21 .  chi  the  other  hand,  if  an  ice  covet 
is  present  and  the  ice  edge  is  receding  downstream  due 
to  melting,  it  will  only  recede  to  the  location  at  which 
the  melting  from  below  equals  the  tendon  \  to  thicken 
from  above,  lhus,  during  periods  of  warming  air  tem¬ 
peratures  when  the  ice  edge  is  receding  downstream, 
the  criterion  of  heat  balance  represented  bv  eg  23  ap¬ 
plies.  During  periods  of  cooling  air  temperatures,  the 
O' l'  isotherm  location  governs  first  ice  appearance  and 
eg  21  applies.  In  the  simulation  that  follows,  the  order 
ot  logic  in  the  calculation  ot  ice  thickening  results  in 
either  the  0"0  criterion  or  the  heat  balance  criterion 
being  naturally  satisfied,  although  the  latter  is  onlv 
approached  asv  mptoticallv  in  time  because  ot  the  un¬ 
steady  nature  of  the  calculations. 
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Figure  4.  Comparison  of  distances  downstream 
to  edge  of  ice  cover  for  the  criteria  of  equilibrium 
heat  fluxes  and  0°C  isotherm  location  as  func¬ 
tion  of  the  air  temperature  and  initial  water 
temperature. 


since  increased  depth  results  in  slower  attenuation  of 
rw  due  to  O  (through  the  velocity  effect)  than  occurs 
by  simple  reduction  of  velocity  as  a  result  of  width 
increase  (since  U  has  small  effect  on  the  attenuation  of 
7"w  ).  The  above  reasoning  is  perhaps  made  more  clear 
by  examination  of  the  decrease  in  Tw  with  distance 
described  by  eq  14,  which  is  presented  in  Figure  2. 
Again,  the  analysis  about  to  be  described  would  allow 
quantitative  examination  of  specific  cases. 

Numerical  simulation 

The  numerical  simulation  is  quite  straight  I  or  ward. 
The  temperature  evolution  along  a  reach  is  computed 
using  eq  5  in  the  difference  form 

7-;+,  =  r;--£— a/  (25) 

w'  •  pcpd 

where  the  superscript  i  denotes  a  time  step,  and  the  sub¬ 
script  j  denotes  a  distance  node.  Since  cq  25  is  in  the 
Lagrangian  sense,  this  requires  that  AXj  =  UfAt.  This 
requires,  in  turn,  that  the  total  reach  be  subdivided  into 
subreaches  with  lengths  varying  if  U  varies  but  this  pre¬ 
sents  no  serious  problem.  Tv pical  time  steps  used  were 
of  the  order  of  1200  to  1800  seconds.  At  the  conclu¬ 
sion  of  each  time  step,  the  thickness  of  the  ice  was 
determined  at  each  /  node  point  using  eq  20  written 
in  the  form 


Figure  3  serves  mainly  to  indicate  the  conditions 
for  which  the  0°C  isotherm  criterion  is  a  poor  approx¬ 
imation  for  the  location  of  the  ice  edge  under  steady 
state  conditions.  The  numerical  simulation  presented 
below  is  sufficiently  simple  computationally  that  the 
unsteady  effects  of  F1  may  also  be  included.  The 
regions  for  which  the  ratio  is  negative,  of  course,  are 
where  the  analysis  presented  herein  is  most  useful,  as 
the  0°C  criterion  would  predict  finite  lengths  of  open 
water  while,  in  fact,  no  open  water  will  exist  except  in 
the  very  near  field.  This  region  also  corresponds  to  the 
conditions  which  Silberman  (1974),  in  a  discussion  of 
Paily  ct  al.’s  (1974)  analysis,  suggested  resulted  in 
plunging  of  the  warm  discharge  beneath  a  cooler,  but 
lighter,  layer  of  near  0°C  water.  Such  a  plunging  phe¬ 
nomenon  undoubtedly  occurs  for  small  velocities  and 
small  densimctric  Froude  numbers  and  clearly  is  the 
case  for  warm  water  flowing  into  large  bodies  of  water. 
Similarly,  the  reappearance  of  open  water  downstream 
from  an  ice-covered  reach  may  be  explained  by  the 
existence  of  a  larger  velocity  downstream  from  a  region 
of  smaller  velocity,  since  0wj  is  increased  more  or  less 
proportionally  to  the  velocity.  The  effect  would  be 
more  likely  to  occur  where  the  lower  velocity  results 
from  increased  depth  rather  than  from  increased  width, 


Ihe  initial  temperature  at  the  upstream  end  of  the 
subreach  was  specified  for  each  time  as  was  the  air  tem¬ 
perature.  For  the  examples  presented,  daily  averages 
were  used  although  it  would  be  quite  easy  to  specify 
a  more  detailed  time  variation.  In  the  calculation  of 
0  at  each  /  node  point,  it  was  necessary  to  use  cither 
cq  6  or  cq  8,  depending  on  whether  the  surface  was 
open  or  ice-covered.  In  the  thickness  calculations,  0; 
was  taken  equal  to  zero  if  Tt  >  0°C.  After  the  new 
ice  thickness  was  calculated,  it  was  set  equal  to  zero  if  a 
negative  value  resulted.  In  the  example  simulations 
presented  here,  /?WJ  and  hu  were  assumed  equal. 

The  problem  also  requires  specifications  of  initial 
conditions.  In  Figure  5  are  presented  two  examples 
of  the  effect  of  specifying  a  constant-thickness  initial 
ice  cover  of  either  zero  or  finite  thickness.  Ta  was 
maintained  at  a  constant  value  of  -5°C.  The  distance 
to  ice  edge  stabilizes,  after  many  days,  for  the  finite 
initial  thickness  of  0.3  m  and,  after  several  days,  for 
the  zero  initial  thickness.  In  the  latter  case,  the  time 
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The  program  is  listed  in  Appendix  A. 


Figure  5.  Example  simulation  of  the  movement 
of  the  ice  edge  in  response  to  an  abrupt  release 
of  thermal  effluent  for  different  assumed  initial 
ice  conditions. 

to  stabilize  is  longer  than  the  time  of  travel  to  the 
equilibrium  ice  edge  location  because  ice  is  grown  until 
the  temperature  front  reaches  the  ice  edge  and  this  ice, 
in  turn,  must  be  melted  before  equilibrium  is  achieved. 

A  better  set  of  initial  conditions  could  probably  be 
obtained  by  running  the  simulation  for  a  few  days  prior 
to  the  period  under  investigation.  Figure  5  also  il¬ 
lustrates  the  potential  errors  in  neglecting  the  eneigy 
required  to  melt  the  ice  cover  in  the  steady-state  analyses 
referred  to  eat  tier. 

There  is  one  particular  circumstance  where  the  above 
analysis  must  be  modified.  When  the  air  temperature 
is  cooling,  the  ice  edge  location  moves  upstream.  How¬ 
ever,  eq  26  is  not  applicable  here  until  an  ice  cover 
exists  and  ice  will  not  form  until  the  water  has  cooled 
to  0°C.  In  this  case,  the  0°C  isotherm  location  is  the 
more  correct  location  than  the  location  given  by  the 
case  corresponding  to  0W,  =  0WJ  as  rj  -*  0  (cq  23  in 
the  steady-state  case).  There  arc  thus  two  locations  for 
the  ice  edge,  depending  on  the  prior  history  of  the  ice 
cover  extent. 

Finally,  while  no  complete  field  data  arc  available  at 
the  time  of  writing  the  simulation  has  been  applied  to 
the  1965  data  presented  by  Dingman  etal.  (1967)  but 
with  an  estimate  of  0,  U  and  D.  The  results  are 
presented  in  Figure  6  and  show  reasonable  agreement 
with  the  three  observations  available. 


Uncertainties  and  limitations 

There  are  a  number  of  uncertainties  in  the  analysis 
presented  above  and  a  number  of  limitations  to  its 
applicability.  Uncertainties,  aside  from  the  difficulties 
of  accurately  representing  the  meteorological  variables, 
include  some  uncertainty  in  the  calculation  of  /twa, 

6wj  and  ha.  For  example,  hwj  is  known  to  increase  by 
about  50%  due  to  the  relief  features  which  form  on  the 
underside  of  the  ice.  Since  these  relief  features  appear 
as  a  consequence  of  <pwj,  it  is  probably  reasonable  to 
use  the  higher  value  in  calculations.  There  is  also  con¬ 
siderable  uncertainty  that  the  ice  edge  location  pre¬ 
dicted  by  use  of  eq  1 9,  which  leads  to  eq  23,  is  appro¬ 
priate.  For  the  case  of  increasing  air  temperatures  and 
a  downstream  receding  ice  front,  it  is  probably  a  good 
approximation.  For  the  case  of  decreasing  air  tempera¬ 
tures  and  an  ice  front  moving  upstream,  it  is  probably 
poorer  because  the  hydrodynamic  forces  acting  on  a 
newly  formed  thin  ice  cover  result  in  initial  thickening 
by  accumulation  rather  than  by  thermal  thickening. 

One  possible  improvement  would  be  to  calculate  the 
ice  production  under  such  cases  and  locate  the  ice 
front  by  use  of  the  equilibrium  accumulation  thickness 
predicted  by  the  analysis  of  Pariset  and  Hausser  (1961 ). 
Even  that  approach,  however,  would  be  limited  by  the 
fact  that  above  a  velocity  of  about  0.6  m  s'1  it  is  dif¬ 
ficult  for  a  cover  to  even  form  by  such  accumulation. 

The  limitations  imposed  by  the  simplifications  of 
the  analysis  are  several.  Longitudinal  mixing  has  been 
neglected  but  this  is  considered  to  have  a  negligible  ef¬ 
fect  compared  with  the  thermal  inertia  of  the  ice  cover 
itself.  The  assumption  of  complete  vertical  mixing 
makes  the  analysis  inappropriate  for  locations  very  close 
to  tl  e  source  of  thermal  effluent,  unless  the  effluent  is 
already  fully  mixed  as  is  the  case  for  reservoir  discharges. 
The  assumption  of  complete  mixing  across  the  width  of 
the  flow,  and  consequent  neglect  of  lateral  mixing  and 
dispersion,  is  important  both  near  the  source  of  the 
effluent  and  well  downstream  near  the  location  of  the 
ice  front.  The  source  problem  is  reasonably  self-evident. 
At  the  ice  front,  the  limitation  occurs  because  of  varia¬ 
tions  in  depth  and  velocity  across  the  width  of  the  flow 
which,  because  of  the  resulting  lateral  variations  in 
dTw/dt,  cause  a  laterally  uniform  temperature  distribu¬ 
tion  to  become  non-uniform.  Thus,  it  is  common  that 
the  ice  front  is  observed  to  be  somewhat  V-shaped, 
with  the  apex  of  the  V  downstream  and  near  the  thal¬ 
weg  of  the  channel.  These  latter  effects  are  currently 
under  analysis  by  introducing  lateral  mixing  into  the 
formulation  and  require  the  additional  consideration  that 
the  dispersion  coefficient  for  closed  surface  flow  is  ap¬ 
proximately  half  that  for  open  surface  flow  at  the  same 
depth. 
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/  njnre  f>.  (  om/sarison  ot  simulation  results  to  ofr- 
st •tuitions  tor  Riverside  lomt ion  in  lop 

figures  are  the  ait  temperature  and  release  i safer 
temperature  variations.  Hot  tom  injure  is  the  simu¬ 
lation  result  I solid  line)  and  ohsen  atimis  ( plotted 
points). 

I  iiully,  the  hydraulics  ut  the  Mow  have  been  as¬ 
sumed  to  tv  steads .  this  is  a  serious  limitation  hut 
could,  in  principle,  he  overcome  by  frequent  use  ot 
numerical  models  ot  open  channel  Mow  to  calculate  the 
applicable  velocities  and  depths  at  vaiious  stages  ot  the 
simulation,  these  models  should  include  the  ettect  ot 
ice  covet  i<n  the  Mow  but  ate  substantially  uncoupled 
hom  the  thermal  analysis. 
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PART  II.  EFFECT  OF  TRANSVERSE  MIXING 
ON  ICE  SUPPRESSION 

Introduction 

Rivers  are  commonly  used  for  the  disposal  of  ther¬ 
mal  wastes  and,  in  most  cases,  the  disposal  occurs  in 
the  form  of  a  side  channel  discharge  of  heated  effluent. 
During  periods  of  ice  cover,  the  effect  of  such  effluents 
is  to  suppress  the  ice  cover  from  its  otherwise  natural 
thickening  and,  in  the  vicinity  of  the  release,  open  water 
often  results.  The  intent  of  this  work  is  to  explore  the 
extent  of  this  ice  suppression  as  a  function  of  the  para¬ 
meters  which  have  most  significant  effect  on  the  sup¬ 
pression.  Part  I  of  this  report  dealt  with  the  ice  sup¬ 
pression  resulting  downstream  from  an  effluent  fully 
mixed  with  the  flow  as  would  be  experienced  down¬ 
stream  from  a  reservoir  release  or  from  a  power  plant 
with  an  outfall  diffuser  operated  so  as  to  provide  com¬ 
plete  mixing. 


Analysis  of  dispersion  and  heat  loss 
The  governing  differential  equation  is 


U) 


where  p  is  the  density  (kg  m‘3),  Cp  is  the  specific  heat 
(/  kg'^C'1),  Tw  is  the  water  temperature  (°C),  Ez  is 
a  transverse  dispersion  coefficient  (m2  s'1 ),  0  is  the 
heat  flux  at  the  top  surface  (W  m‘2),  D  is  the  depth  (m), 
x  is  the  longitudinal  coordinate,  z  is  the  transverse  co¬ 
ordinate,  U  is  the  mean  velocity,  and  t  is  time.  Equation 
1  implies  that  mixing  of  the  effluent  with  the  flow  is 


complete  over  the  depth  and  that  there  is  no  lateral 
mixing  due  to  transverse  velocities  (secondary  currents). 
The  effect  of  longitudinal  dispersion  is  also  neglected, 
although  for  rivers  it  is  negligible  compared  with  the 
transverse  mixing  effects  and  the  unsteady  effects  due 
to  variation  in  0  resulting  from  air  temperature  varia¬ 
tions  with  time.  The  product  pCp  is  effectively  con¬ 
stant  over  the  range  of  7"w  of  interest  so  that  eq  1  may 
be  written  in  the  form 


The  heat  flux  0  at  the  surface  depends  on  whether 
or  not  an  ice  cover  is  present.  If  an  ice  cover  is  present, 
0may  be  reasonably  represented  by  an  expression  of 
the  form 

^wi  ~^wi  (3) 

where  0w|  is  the  heat  flux  from  the  water  to  the  ice, 

/twj  is  a  heat  transfer  coefficient  applied  to  the  tempera¬ 
ture  difference  rw-7"m  between  the  flow  and  the  ice/ 
water  surface  which  is  at  a  temperature  7"m  =  0°C  be¬ 
cause  of  the  state  relation.  The  heat  transfer  coefficient 
hwt  depends  on  the  flow  variables.  By  analogy  with 
closed  conduit  turbulent  heat  transfer,  hwj  is  deter¬ 
mined  from  a  Nusselt-Reynolds-Prandtl  number  cor¬ 
relation  of  the  form 


where  />wj  R/kw  is  the  Nussclt  number,  URp/p  is  the 


Reynolds  number,  pCp/Aw  is  the  Piandtl  number,  R 
is  the  hydraulic  radius,  A'w  is  the  water  thermal  conduc¬ 
tivity  (W  m  1  °C-1 ),  and  p  is  the  dynamic  viscosity 
s'  ).  C  is  somewhat  uncertain  but  is  of  the 
order  of  0.023  to  0.030  depending  upon  the  under¬ 
surface  roughness  of  the  ice.  I  he  Prandtl  number  for 
water  is  13.6  at  0  C  and  decreases  with  increasing  tem- 
Pfature.  At  0°C,  dPr/dlw  =*  -0.4  “c'.so  in  the  range 
0  C  to  4°C  the  error  in  assuming  Pr  constant  yields  at 
most  about  5%  error  in  calculation  of  /?w|  and  will  be 
neglected  in  view  of  the  uncertainty  in  C. 

If  the  water  surface  is  open  to  the  atmosphere,  then 
0  depends  on  the  energy  budget  at  the  water  surface, 
for  present  purposes,  we  will  approximate  0  by 


the  flow  velocity  is  unclear.  In  open  water  areas  far 
from  the  ice  edge,  it  is  reasonable  to  expect  E  l 
while  under  the  ice  cover  far  from  the  edge  L  ,  l  , 

I  he  nature  of  the  numerical  simulation  presented  below 
makes  it  convenient  to  assume  an  abrupt  change  in  / 
at  the  ice  edge.  I  he  finite  difference  algorithm  used 
effectively  causes  a  transition  from  L  to  to  /  over  a 
distance  depending  on  the  width  of  the  grid  elements. 

Analysis  of  ice  thickening  and  melting 

I  fie  melting  and  thickening  of  the  ice  cover  aie 
governed  by  the  energy  balanceat  the  water  ice  .n'.ei 
face  (I  ig.  I ) 


^W4  “  ^  W  4) 


where  0WJ  is  the  heat  flux  from  the  water  to  the  air, 

is  a  heat  transfer  coefficient  applied  to  the  difference 
between  the  water  temperature  Tw  and  the  ambient  air 
temperature  Ta.  Coefficient /?WJ  has  large  dependence 
on  the  wind  velocity,  generally  of  the  form 

hWi  =  <?+/>  I/W  (6) 

where  a  is  the  heat  transfer  coefficient  for  still  air  and 
b  V*  is  the  wind  effect  with  l'w  the  ambient  wind 
velocity  at  some  specified  distance  above  the  surface. 

In  numerical  examples  presented  later,  /?WJ  will  be 
taken  constant  although  the  effect  of  eq  6  may  easily 
be  included  in  the  numerical  analysis  (or,  for  that  mat- 
ter,  /?WJ  can  be  calculated  from  even  more  detailed 
energy  budget  calculations  of  0W4). 

The  transverse  mixing  coefficient  also  depends  upon 
whether  or  not  the  flow  is  ice-covered  or  open.  We  here 
relate  Et  to  a  coefficient  k  times  the  product  of  the 
shear  velocity  U .  and  hydraulic  radius  R  in  the  form 

E,=kU.R  (7) 

where  for  open  channel  flow  R  =*  0  and  for  covered 
flow  R  =»  0/2.  Denoting  Lt  for  open  flow  by  Elo  a  -d 
for  covered  flow  by  E jj(  then 


where  </>.  is  the  heat  flux  by  conduction  through  the  i«.e 
cover  (W  m  ’),  p,  is  the  density  ol  ice  (kg  m'J),  X  is  the 
heat  of  fusion  (/  kg  1 ),  q  is  the  ice  thickness  (m).  and 
di}/dt  is  the  late  of  thickening.  Ihe  conductive  tlu\ 

0j  through  the  ice  cover  is  treated  in  a  quasi-steady 
manner  by  assuming  a  linear  temperature  profile  thiough 
the  ice  thickness;  thus, 


where  k{  is  the  thermal  conductivity  of  the  ice  (W  m  1 
(  ),  and  /  4  is  the  top  surface  temperature  ot  the  ice 

cover  and  we  require  that  <  0°C  because  ot  the  stale 
relationship.  In  the  absence  of  evaporation  01  conden- 
sation  on  the  top  surface,  0,  =  <fiu,  where  0U  is  the  heat 
flux  from  the  surface  to  the  atmosphere.  0U  may  be 
calculated  in  a  manner  similar  to  0WJ  thiough  intio- 
duction  of  a  heat  transfer  coefficient  hu  applied  to 
the  difference  7\- Tt  in  the  form 

(Ti-Tt).  (II) 

Using  0 j  0U  allows  us  to  eliminate  1\  between  eq  10 
and  1 1  and  results  in 


Elo  -2Eti  k  0.0. 


Engmann  (1977)  has  measured  Et  for  both  covered  and 
open  flows  and  found  results  consistent  with  eq  8.  with 
k  in  the  range  0. 1 5  to  0.2  which  may  be  compared  to 
a  commonly  used  value  of  0.23  (Okoye  1970).  How¬ 
ever,  if  the  river  is  not  straight,  considerably  larger 
values  of  k  are  observed  (Paily  and  Sayre  1978). 

The  situation  for  a  partially  open  ice  cover  where 
the  edge  of  the  ice  cover  is  more  or  less  aligned  with 


0.=I=di 


We  could  also  in  a  similar  manner  add  the  effects  ol  a 
snow  cover  by  introducing  an  additional  term  » it/ks  in 
the  denominator  of  eq  1 2. 

Substitution  of  eq  1 2  in  eq  9  then  results  in 
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Figure  I.  Definition  sketch  for  heat  transfer  analysis. 
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which,  if  Tt  and  £w  arc  constant,  may  be  readily  inte¬ 
grated.  They  seldom  are,  however,  so  in  the  numerical 
analysis  eq  13  will  be  integrated  in  step-wise  fashion 
using 


Ar,=  ^L 
p/A 


T"-T±-h 


jl+_L 


wi 


(K-TJ 


04) 


-__±_ 

PCPD 


(16) 


where  the  £  values  depend  on  whether  or  not  the  flow 
is  ice-covered  (£  =  £I()  or  open  (£  =  £io)  at  the  respec¬ 
tive  /  grid  points.  The  <p  term  is  also  dependent  on 
whether  or  not  the  flow  is  ice-covered  and  we  choose 
the  /',/'  location  for  the  determination.  Thus 


0wj  =  hwi  (ice  cover  P^cnt)  (1 7) 

0W1  =£Wi  (Tj-T,)  (open  surface).  08) 


and,  of  course,  if  the  resulting  value  of  n  <  0  after  a  Solving  eq  16  explicitly  for  7f '  results  in 

time  step,  then  rj  is  set  equal  to  0  and  open  water  exists. 


Numerical  simulation 

Because  the  numerical  simulation  presented  below 
considers  the  effect  of  variations  in  the  initial  effluent 
temperature  with  time,  it  is  convenient  to  transform 
eq  2  to  a  Lagrangian  coordinate  system  moving  at  the 
mean  velocity  of  the  flow.  Thus,  introducing  Udt  =  dx, 
eq  2  becomes,  for  a  parcel  of  water  moving  with  the 
flow, 


7j+ 1  =  tf  t{_,  *b  rj+c  7"j+ 1  +cf  7",  (19) 

where  the  coefficients  a,  b,  c  and  d  are  given  by 

—  -(£j*£M>  (20) 


6=1- 


2(A/)2 

Af6wi  fa 


PCp°  2(Az)2 
(ice  cover) 


(£M+2£i+£,  ,) 


(21) 


The  explicit  finite  difference  equation  used  in  the  simu¬ 
lation  is 


Af 

=  (fi*i*fi>  ri'*i-<fi*i*2fi*fi-,>  ^(fi^i-i)  y|-» 
2  (A/)2 


6  = 


1-  ^2L-_£L_(£  , 

f>CP°  2(A/)2  ' 

(open  surface) 


+2^£h) 


c-_^(£i+,+£|) 
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Figure  2.  Results  of  example  simulation  for  a  smaller  river  under  steady-state 
conditions. 
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Figure  3.  Lateral  variation  of  water  tempera¬ 
tures  for  example  simulation  after  five  days. 


</  =  0 

(ice  cover) 

(24) 

pCpD 

(open  surface). 

(25) 

Again,  the  value  of  the  t 's  depends  on  whether  or  not 
an  ice  cover  or  an  open  surface  is  present. 

The  above  explicit  scheme  is  stable  for 

2fAr<, 

(A/)2 

(26) 

and  the  appropriate  value  for  £  is,  of  course,  E0, 
it  is  larger  than  £(  for  the  same  total  depth. 

since 

Initial  conditions  must  also  be  specified.  Since  one 
seldom  has  detailed  knowledge  of  ice  thickness  varia¬ 
tions,  the  choice  is  practically  between  an  initial  uni 
form  ice  thickness  or  an  initial  open  water  condition. 

In  either  case,  even  for  steady  effluent  temperatures 
and  steady  air  temperatures,  it  takes  some  time  for  an 
equilibrium  state  to  be  reached.  In  the  case  of  an  initial 
zero  thickness  ice  cover,  this  is  due  to  the  requirement 
of  melting  ice  which  has  grown  during  the  period  before 
the  thermal  “front”  reaches  a  downstream  point.  In 
the  case  of  an  initial  uniform  ice  thickness,  the  time 
to  reach  equilibrium  may  be  quite  long  because  of  the 
large  thermal  inertia  of  the  ice  cover.  In  the  example 
simulations  which  follow  the  initial  conditions  have 
been  taken  as  zero  thickness  and  -  0°C.  In  applica¬ 
tion,  the  most  reasonable  approach  would  seen)  to  be 
to  run  the  simulation  for  several  days  prior  to  the  period 
of  interest  unless,  of  course,  the  situation  to  be  in¬ 
vestigated  is  the  one  of  abrupt  release  of  a  thermal  ef¬ 
fluent. 

A  complete  listing  of  the  computer  program  may  be 
found  in  Appendix  B. 

Example  simulations 

In  this  section  two  example  simulations  are  presented. 
The  first  example  is  that  of  a  smaller  river  50  nr  wide 
and  2  m  deep  with  a  mean  velocity  of  0.5  ms’.  1  he 
thermal  effluent  is  distributed  over  an  initial  width  of 
10  m  at  a  temperature  of  1.0'’C,  which  corresponds  to 
a  heat  load  of  42  MW.  An  initial  thickness  of  0.0  m  is 
assumed,  and  the  air  temperature  is  maintained  con¬ 
stant  at  -5°C.  Other  paianreter  values  used  are 
=  25  W  m2°C'’  (corresponding  approximately  to  a 
mean  wind  velocity  of  about  4.5  ms1),*  0.2,  U* 

=  0.04  m  s'’ ,  and  C  =  0.023.  A  grid  spacing  of  A z 
=  10  m  and  Ax  300  m  (Af  =  bOO  s)  was  used. 

In  I  igure  2  the  ice  edge  location  is  shown  for  times 
of  one,  two,  three,  four  and  five  days  after  initial  ther¬ 
mal  effluent  release,  and  it  is  seen  in  the  figure  that 
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Figure  5.  A  typical  winter  daily  air  temperature  wriation  used  in  the 
example  simulation  of  a  wide  river. 


equilibrium  was  achieved  after  about  three  days.  Also 
shown  in  Figure  2  is  the  boundary  of  the  region  unaf¬ 
fected  by  the  thermal  effluent  after  five  days  (defined 
by  an  ice  thickness  differential  of  Arj  <  0.001  m  for 
»j  =  0.1 1  m)  and  the  location  at  which  the  ice  thickness 
was  one-half  its  undisturbed  value  after  five  days.  In 
Figures  3  and  4  are  presented  lateral  profiles  of  water 
temperature  fw  and  ice  thickness  r?  respectively,  at 
various  distances  downstream  from  the  release  point,  in 
both  cases  for  t  =  5  days.  We  also  note  that  had  the 
thermal  effluent  been  fully  mixed  across  the  width  of 


the  river  at  *  =  0  the  ice  edge  would  have  been  located 
in  the  vicinity  of  1.5  km. 

The  second  simulation  is  for  a  very  wide  river,  in 
fact  sufficiently  wide  that  reflection  from  the  far  hank 
has  no  effect  on  the  results  both  because  of  its  distance 
and  the  decrease  in  Fw  due  to  heat  losses  to  the  ice 
cover.  The  flow  was  characterized  by  a  depth  of  4  m 
and  a  mean  velocity  of  0.5  ms-',  and  a  heat  load  of 
540  MW  was  distributed  over  an  initial  width  of  32  m 
with  a  corresponding  water  temperature  of  4.0°C.  lhis 
heat  load  was  maintained  constant  while  the  air  tem- 
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Figure  6.  Variation  of  length  of  open  water 
for  the  wide  river  e sample  simulation. 
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Figure  8.  Maximum  and  minimum  open  water 
extents  for  the  wide  river  example  simulation. 
perature  varied  as  shown  in  Figure  5  (taken  from  an 
actual  temperature  record  in  the  mid  western  United 
States  for  1972).  An  initial  thickness  of  0.0  m  was 
used  and  othe*  parameters  were  as  follows:  hWJ  =  25 
W  m  J  °C  \  k  =  0.2,  U .  =  0.04  m  s  ',  C  =  0.03,  A/ 

=  16  m,  Ax  =  1 350  m  (corresponding  to  Af  =  2700  s). 

Figure  6  presents  the  resulting  length  of  open  water 
and  Figure  7  the  maximum  width  of  open  water,  the 
location  of  which  varied  somewhat  in  the  longitudinal 
direction  (from  about  6  km  to  25  km).  The  extreme 
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Figure  7.  Variation  of  maximum  width  of  open 
water  for  the  wide  river  example  simulation. 

cases  of  open  water  extent  corresponding  to  24  lanuary 
(maximum  extent )  and  28  lanuary  (near  minimum 
extent)  are  presented  in  Figure  8.  The  extreme  fore¬ 
shortening  of  the  x  scale  relative  to  the  z  scale  in  Figure 
8  (a  ratio  of  250:1 )  somewhat  obscures  how  long  and 
narrow  are  the  open  water  extents.  Finally  we  note  that 
the  example  presented  is  also  valid  for  a  release  at  the 
centerline  with  a  heat  load  of  twice  that  used  since  the 
boundary  condition  at  the  bank  (dT/dz  =  0)  is  the 
same  by  symmetry. 

Field  comparison 

It  is  desirable  to  test  the  simulation  against  actual 
field  observation*.  One  very  limited  set  of  data  was 
obtained  of  the  open  water  that  existed  downstream  of 
the  Riverside  Power  Plant  on  the  Mississippi  River  up¬ 
stream  of  Bettendorf,  Iowa.  The  data  consisted  of  ob¬ 
lique  angle  aerial  photography  of  the  open  water  area 
and  were  obtained  on  17  February  1979.  Figure  9 
shows  a  view  of  the  power  plant  and  near-field  mixing 
/one.  On  the  cover  there  is  a  view  looking  downstream 
from  the  vicinity  of  the  plant.  Figures  10  and  1 1  are 
similar  views  from  points  farther  downstream.  In  Figure 
1 1  an  ice  edge  may  be  seen,  although  farther  down¬ 
stream  there  were  numerous  patches  of  open  water  in 
the  plume  area.  Other  data  were  obtained  subsequent 
to  the  aerial  photography;  their  sources  are  as  follows: 

River  hydrography .  5  m  depth  fairly  uniform 
along  west  side  of  river;  obtained  from  hydro- 
graphic  charts  of  the  river. 

Mean  velocity.  0.61  ms'1,  obtained  from  a 
float  measurement  on  13  February  1979.  River 
flow  discharges  were  reasonably  uniform  over  the 
period  13  to  17  February. 
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f  mart-  10.  I  lew  t  >f  open  water  \e\eral  kilometers  downstream  from  the  Kherside  power  plant. 


Figure  11.  View  downstream  near  end  of  open  water  extent.  (Our  apology  for  the 
poor  quality  of  the  photograph.) 


Daily  air  temperatures',  average  of  daily  maxi 
mum  and  daily  minimum  air  temperatures  re¬ 
corded  at  the  Quad-City  Airport  in  Moline,  Il¬ 
linois,  which  is  located  about  six  miles  from  the 
river  site. 

Initial  width  of  mixing-,  determined  from 
aerial  photography  to  be  approximately  30  m. 

Effluent  source  temperature.  2.97°C.  This 
value  was  obtained  using  design  temperature 
rise  of  water  passing  through  the  various  plant 
units  after  correcting  for  plant  capacity  (figures 


supplied  by  lowa-lllinois  Gas  and  Electric  Com¬ 
pany).  The  effluent  water  was  then  mixed  over 
a  width  of  30  m  of  the  river  flow  to  obtain  the 
2.97°C.  Considerable  uncertainty  exists  in  this 
value  since  only  monthly  plant  operating  figures 
were  available.  Further,  the  intake  and  outfall 
geometry  are  sufficiently  complicated  that  a  sim 
pie  near-field  mixing  analysis  was  not  feasible. 

Wind  speed'.  Not  recorded.  In  lieu  of  site  in¬ 
formation,  /twa  and  hjj  were  arbitrarily  taken  as 
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Figure  1 2.  Comparison  of  numerical  simulation 
results  for  16,  17,  18  February  and  observed  open 
water  extent  on  17  February. 

25  W  m  *  °C‘' ;  this  corresponds  approximately 

to  an  average  wind  speed  of  4.5  ms’. 

Using  the  above  information,  three  simulations  were 
run  on  the  computer  using  values  for  the  mixing  co¬ 
efficient  of  *  =  0.20,  0.40  and  0.80  m2  s'.  The  last 
value  provided  the  best  agreement  with  observation  of 
the  three  runs  and  the  resulting  open  water  extent  is 
shown  in  I  igure  1 2  for  the  days  of  1 6,  17  and  1 8 
February.  Also  shown  is  the  observed  open  water  ex¬ 
tent  on  1 7  February. 

Some  interpretation  of  Figure  12  is  in  order.  Ap¬ 
proximately  7  km  downstream  from  the  plant  the  ice 
edge  appears  to  be  pinched  at  the  location  of  a  sus¬ 
pension  bridge  over  the  river.  This  is  perhaps  explained 
by  the  presence  of  a  pier  located  approximately  at  the 
ice  edge  location  which  may  act  to  stabilize  the  cover . 
At  a  distance  of  II  km  Lock  and  Dam  1 5  severely 
changes  the  local  flow  geometry  since  all  the  river  flow 
passes  through  roller  gates  on  the  opposite  side  of  the 
river.  Evidence  of  the  thermal  plume  existed  down¬ 
stream  of  L  &  D  1 5  but  no  attempts  were  made  to 
simulate  the  complex  flow  conditions  resulting  at  the 
dam.  The  simulation  was  arbitrarily  extended  using 
the  upstream  depths.  The  "ice  edge"  shown  in  I  igure 


12  is  that  also  depicted  in  Figure  11.  However,  ex¬ 
tensive  patches  of  open  water  existed  in  the  thermal 
plume  area  downstream  and  the  solid  line  shown  in 
Figure  I  2  is  the  outline  of  this  area.  The  accuracy  of 
the  simulation  is  difficult  to  evaluate  from  these  limited 
data,  both  because  of  the  uncertainty  in  the  input 
parameters,  and  because  ot  the  rapid  advance  upstream 
of  the  ice  edge,  resulting  from  a  four-day  cold  period. 

The  simulation  was  performed  using  a  transverse  grid 
interval  of  I  5  m,  a  total  grid  width  ot  1  50  rn,  and  a 
time  step  of  400.8  s. 

It  is  planned  to  conduct  a  more  detailed  field  study 
at  this  site  to  include  direct  measurement  ol  water  tern 
perature  profiles,  ice  extent  and  meteorological  variables. 
These  data  will  allow  a  better  determination  of  the 
quality  of  performance  of  the  simulation. 

Uncertainties  and  limitations 

The  uncertainties  present  in  the  analysis  largely  re¬ 
late  to  the  magnitude  of  the  coefficients  used.  and 
hu  (assumed  equal  in  this  analysis)  can  always  be  better 
approximated  if  more  extensive  meteorological  infor¬ 
mation  is  available.  The  magnitttde  of  C  used  to  esti¬ 
mate  bwi  should  be  verified  by  actual  field  measure¬ 
ments  and  one  virtue  of  the  present  work  is  a  means  of 
estimating  whether  a  particular  field  site  is  sufficiently 
well-mixed  horizontally  that  measurements  of  temper¬ 
atures  in  the  streamwise  direction  will  accurately  \  ield 
a  value  of  C.  The  ice  edge  location  results  from  a 
melting  (or  freezing)  analy  sis  that  assumes  the  cover 
is  stable  for  very  small  values  of  r).  This  is  probably 
reasonable  in  the  lateral  direction  since  at  the  ice  edge 
the  thickness  changes  rapidly  with  z.  It  is  less  certain 
at  the  downstream  edge  since  the  ice  thickness  tends 
to  “feather"  and  be  quite  thin  over  tairly  large  x  dis¬ 
tances.  Other  limitations  in  location  of  the  downstream 
extent  are  discussed  in  Part  I. 

The  analysis  has  assumed  a  uniform  distribution  of 
velocity  and  constant  depth  in  the  z  direction  which 
is  particularly  convenient  since  time  steps  then  cor¬ 
respond  to  distance  steps.  Holley  et  al.  (1972)  have 
discussed  the  effect  of  lateral  variations  in  velocity  and 
depth  on  the  dispersion  coefficient,  and  it  can  be  sig¬ 
nificant.  One  improvement  in  the  present  analy  sis 
would  be  to  incorporate  such  effects,  perhaps  even 
with  a  quasi-steady  shift  of  the  velocities  towards  the 
open  areas  as  the  cover  is  melted  away.  One  procedure 
possible  is  that  suggested  by  Say  re  and  t  eh  (1973),  who 
utilized  the  Manning  relationship  on  a  local  basis  to  de¬ 
termine  the  transverse  distribution  ot  water  discharge. 
Finally,  secondary  currents  can  result  in  considerable 
magnification  of  the  effective  value  ot  k  when  the 
river  meanders.  Incorporation  of  these  effects  would 


probably  require  an  implicit  algorithm  rather  than  the 

conveniently  simple  explicit  algorithm  used  herein. 
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APPENDIX  A:  UNSTEADY  FULLY  MIXED 
ICE  SUPPRESSION 
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0  ASHTON  APPtl  1978  UNSTEADY  mit.Y  ‘I*F0  Iff  siiRPRFSSTOr 
r H T S  PROSPAM  COMPUTES  DOWNSTREAM  FVOLUTTOfc  or  AN  IFF  COVER 
in  PfSP.V'SF  TO  PflFASE  Of  4  fltl.1  Y  I'TTED  TwfrvAI  f  f  f  I  "  F  T 
AND  VARYING  AIR  TE*PJ  PTVfo  flOw  TS  •  SSII**f  0  TO  HF  STFADy 
»L#0#U#FTA,  TWOOT  ,  u  ‘RE  I  E“,C.Th,1FPTM,ve  »f.  VELOCITY, 

ICE  THICHMESS.wATFR  T  F  "»PF  R  »  TU»F  ,  tun  WIDTH  FOR  FACH  $  'IM  P  F  ft  0  H 
OAT  Aon  «T‘J  APf  n  A  T I  Y  ATP  TEMPSAMD  FfFIlif.'T  SOlIPCF  TFWPS 
0 THE MS  ton  SALC6C>»SF»(/>O),S0(80>,RAT(80>,STN(/>n) 

DIMENSION  A I  (100)  ,0(100)  ,0(100)  ,ETA  (100)  ,  TWOI'T  flOli)  ,H(  100) 

DI*FNSI0N  CV’AM(70)»V(70) 

C'-lAR  AND  V  APf  DATIY  H  F  A  T  TRANSFER  COFF  8  '.'JMD  (PffS 
PfAl)  900»NT,NSR,Al  TOT, DISC* 

FORMAT (?I10,?E10.0) 

PRINT  Ml/., NT  ,N$R,ALTOT,OISCM 

FOP WAT/1M1 ,Sx, »HT  NSR  Al  TOT  o ( SCh  • ,7T4,7f 1?.8) 

Al  TOT  IS  IENGTM  Of  STUDY  REACH  ( W )  ,  OTSfM  IS  DISCHARGE  0*3/S>, 

NSR  TS  mumper  OF  SOH0E*CmFS  foo  WHICH  V  A I  Ilfs  OF  width  OF 
CHANNEL  A  HO  DFPTH  APf  AVAJLAKLF 
IT  IS  MIIMHFP  Of  DAYS  OF  ST“o(  aTTON 
PEAD  901.(Si»(J).SP<J)*S*L<J)»J»1»0SP) 
format  (JFin.i) 

PRINT  R17#<SR(J),S0<J),SAL(.1>/J«1,NSR> 

FOPWATd*.'  SP  a  '/F10. 3,*  SO  *  ,»F10.3>*  S  A I  a  ',F10.3) 

SO  AND  SD  APE  width  AND  DFPTm  at  SUCCESSIVE  *  DISTANCES  sal 
SUBDIVIDE  TOTAI  PFACH  INTO  SIIHPFAOMFS  I.  T  T  H  Df(T  TPAVFI.  T I  Hp 

dflt  *  3000. 

CALFOl  ATF  ml  AH"  TOTAI  rPAVEI  TIME 
T R  A  V  *  0.0 
no  10  I  ■  1  ,  N  $  0 

TPAV*TPAVASP(J)ASD(J)*SAL(J)/OI«rH 
A'LaTPAV/OFl  T 
TRAVsTRAV/3600. 

PRIST  810.TRAV 

FORMAT (  1  0  X  , 1  TOTAI  TIME  Of  TRAVEL  »  '.Flo.?.'  HOOPS') 

CALCUIATF  MFU  si'ORFATm  |  FfJOTHS  .IT  DFl  T  Tt^-E  Of  toavFI. 

SONS  A 1  *SAL (1  ) 

SllHALaO<0 

I»1 

DO  15  Jal.Nl 

U(J)*DISCM/(SH(T)*SD(I>) 

AL(J)«Df  LTai.KJ)  A  . 

0<J)aSD<I) 

S.IHACSLIMACALO)  THECiOlWG  PaaGE  BUWK  -  NCr  n, am 

0EL8L-SUMAL-SUMSAL  ...  H^JAAYAV  -  mj  i  Jl.U'ifcD 

If <DFL»L>1S,1S,1? 

I»I*1 

Sll**SAL*SOMSAL  ASAl  (I) 

CONTtNOE 

routine  Omit  rood  if  all  sai  i  ongf®  tha»  dflt*m(j) 

PRINT  811 

FORMAT  MX,  •  AI.CJ  >  *,3X,  'tJ(J)  '  ) 

PRINT  813 

FORMAT  M», '  MFTFPS ' ,3X, *«/s ' > 
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A  facsimile  catalog  card  in  Library  of  Congress  MARC  format  is 
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